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SAMPLING THEOREMS ON BOUNDED DOMAINS 



MASSIMO FORNASIER AND LAURA GORI 



Abstract. This paper concerns with iterative schemes for the perfect reconstruction of func- 

■u tions belonging to multiresolution spaces on bounded manifolds from nonuniform sampling. The 

'^l,( I schemes have optimal complexity in the sense that the computational cost to achieve a certain 

fixed accuracy is proportional to the computed quantity. Since the iterations converge uniformly, 

one can produce corresponding iterative integration schemes that allow to recover the integral 

of functions belonging to multiresolution spaces from nonuniform sampling. We present also 

an error analysis and, in particular, we estimate the L^-error which one produces in recovering 

^^ , smooth functions in //*, but not necessarily in any multiresolution space, and their integrals 

I ^ 1 from nonuniform sampling. Several uni- and bi-variate numerical examples are illustrated and 

f^ l ' discussed. We also show that one can construct a rather large variety of multiresolution spaces 

on manifolds from certain refinable bases on the real line formed by so-called GP-functions. 

This class of functions that contains in particular B-splines has remarkable properties in terms 

C^ ■ of producing well-conditioned bases. The resulting multiresolution analyses are well-suited for 

the application of the iterative recovering of functions from nonuniform sampling. 

AMS subject classification: 65D05, 65D15, 65D32, 65T60, 94A20 

^ ' Key Words: Nonuniform sampling, quasi-interpolation, multiresolution analysis, banded matri- 

^^ ' ces, iterative integration formulas. 

(N' 

^ ' 1. Introduction 

^. . 

f^ , The reconstruction of a signal from sparse and nonuniform sampling data is a well known 

^O ' problem |4, and mainly addressed to guessing or learning some relatively small missing part of 

a function by using the information of the relevant known part. The problem of learning from 
r^ ' examples is a corner stone in information theory and artificial intelligence. Very recently a prob- 

jrt , ability theory of what can be learned from a relatively large distribution of data according to 

Cj ' a (unknown) probability measure has been formalized in the beautiful contribution by Cucker 

.. I and Smale 10 . Inspired by this work, subsequent papers, for example jS|, illustrate methods to 

^ ' construct estimators of the probabilistically best solutions of the learning problem. 

Jv>( , Besides the statistical learning, mathematical methods and numerical algorithms based on de- 

5^ ' terministic techniques have been developed to compute missing parts of signals from few and 

Cd I sparse known sampling information and we refer, for example, to classical works of Feichtinger et 

al. ^^ElEin and Grochenig et al. |21|^. These methods and algorithms are essentially based 
on a quasi-interpolation of the function defined on the Euclidean space by means of suitable series 
expansions of irregularly shifted (translated) basic functions or, in its discrete version ^I^ IH], by 
discrete finite series expansion of complex exponentials. One of the typical applications is of course 
the restoration of digital signals and images, where some samples (ndr. pixels) can be corrupted, 
e.g., by scratches as it happens in old movies, or even completely missing, as in case of time-series 
of astronomical observation data. As a nice example of the use of such techniques for real appli- 
cations, in ,yi one of the authors discusses nonuniform sampling methods in combination with 
suitable variational models in order to restore colors of missing parts of destroyed art frescoes from 
those of few known fragments (previously detected and placed in their original sites |18II19| ). and 
the gray levels of the missing parts due to some pictures of the frescoes taken prior to the damage. 

In this paper we want to present sampling theorems for functions defined on compact subsets 
Q, of M''. Because of the boundedness nature of the domain such theorems are nicely suited for 
applications to signals in concrete situations. As we will show, even the proofs result relatively 

1 



2 MASSIMO FORNASIER AND LAURA GORI 

elementary with respect to those for functions defined on the whole Euclidean space ^1 |2]- 
The reconstruction from nonuniform sampling of functions defined on intervals has been also 
considered in |15[ I24| . In JS] the functions are modeled as complex trigonometric polynomials 
and in |21] as restrictions of functions defined on R as linear combination of shifted compactly 
supported functions. None of these two approaches can be really extended to the reconstruction 
of functions defined on manifolds, e.g., on the sphere. The approach that we present in this 
paper can be straightforward generalized to compact manifolds with or without boundary and it 
is computationally very efficient as the methods in |15[ I24| . In the following, the functions are 
assumed to belong to level spaces of a multiresolution analysis directly constructed on the domain 
(or manifold) [HIIZIEIEIj and not as a restriction of functions defined on larger spaces. Moreover, 
the fact that one can associate wavelet bases to these multiresolution analyses to characterize 
function spaces of smooth functions, e.g., Sobolev spaces, allows to derive certain error estimations 
that would not be available otherwise. 

The paper is organized as follows: Section 2 describes the abstract setting and the general 
requirements of the multiresolution spaces we use to model functions, and the corresponding 
sampling theorems for the reconstruction of functions from nonuniform sampling. Section 3 is 
devoted to the numerical implementation of the iterative reconstruction formulas derived in Section 
2. Several numerical examples are shown and discussed. In Section 4 we show how the iterative 
formula can be also used to implement iterative numerical integration schemes from nonuniform 
sampling. We derive error estimators for smooth functions belonging to certain Sobolev spaces in 
Section 5. The construction on compact manifolds of multiresolution analyses with the properties 
that allow the application of the sampling theorems is addressed in Section 6. We will show that 
one can construct a rather large variety of multiresolution spaces from certain refinable functions 
on the real line chosen in the class described in ,22,. 
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2. Elementary axiomatic sampling theorems in LP{i}) 

Let us assume that 
MRAl) is given a multiresolution analysis of finite dimensional spaces V = {V^Ijgn such that for 



all p G (1, oo) 



i''(o) 



xneVo{n)cVi{n)c...v,{n)cLP{n), UT/,(f^) ^LP{n); 

j 

MRA2) each space has a suitable basis of nonnegative continuously differentiable functions $j :— 
{(j)j^o,- ■ ■ ,(t>j,Nj} C C^(ri,M+), such that $j_i = Aj^j where Aj e Mnj^ixn^O^) is a 
suitable (stochastic) scaling matrix, '2^~ J2k=o^3-k = 1 f^'' ^^1 J ^ I^- Moreover, we 
assume that diam(supp((/)j,fc)) < C2^-', #{fc : supp{(j)j^k) H supp((/)j_fe') y^ 0} < m for some 
m G N, uniformly with respect to j, and ||V(/)j.fc||oo < C2^^'^^'^^^\ uniformly with respect 
to both k, j; 

MRA3) Pj : C{n) -^ Vj{n) is a bounded linear projector, i.e., Pjf ^ f for aU / e Vj{n), for aU 
j £ N. Moreover we assume that 

(1) / - lim P,f, for aU / e LP{n), 

with convergence in U'{^). 
Definition 1. A set X = {xi}^^q of sampling nodes in Vt is called A-dense if 

(2) 



17= l\jBs{x,)\{^n, 
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for < (5 ;= A^^, where Bs{x) is the ball centered at x of radius 6. 
One has the following immediate result. 



Lemma 2.1. For a given A-dense set X of sampling nodes, there exists a set of functions ^ = 
{ipe}eLo ^ L°°{fl) with the following properties 

(1) < ^f < 1 for all £ = 0, ..., M; 

(2) supp(V'^) C Bs{xi); 

(3) T.lLo^i = 1- 

Associated to a A-dense set X of sampling nodes one can define the following quasi-interpolation 
operator given by 

(3) Q*,x/ :- Y. /(^^)^^' for all / e C{9). 

In the following we will make use also of another quasi-interpolation operator. For any j G N and 
for all k — 0, ...,Nj, let us choose a sampling point £^j^k G supp((/)j,fc). Therefore, associated to the 
sampling set Sj = {^j.k}k=o ^^ consider the following operator 

(4) 5$,H,/ := 2-*/2 J2 f{^i.k)4>j.k, for all / G C{n). 

e=k 

Let us formulate now the main result of this paper. 

Theorem 2.2. If V — {VjjjeN is a multiresolution analysis in L^lfl) for p e (l,oo) with prop- 
erties MRAl)-MRA3), then for any j G N there exists Aj > large enough such that for any 
Aj-dense sampling set Xj := X{Aj) any function fj G Vj{n) can be exactly recovered from its 
samples {fixjj)}^Ji by the following iterative algorithm 

(5) /]"+'' = p,Q*„x,(/,-/]"') + /]"\ n>i, /f -p,g*,,x,/,. 

In fact one has 

(6) /, := lim /j"\ 

where the convergence is uniform on fl, and then it is valid in L'^iTl) for all q G [1, cxd]. Moreover, 
for all f G LP(fl), let fj = Pjf, one has 

(7) /= lim lim f^"\ 

j — *oo n — 'oc -^ 

being / defined by the iterative scheme (0 by means of suitable Aj-dense sampling values 

The iterative scheme Q is analogous to that appearing in pQ, but here it has been adapted for 
the reconstruction of functions defined on bounded domains. 

REMARK: For / G C(0) but / ^ Vj{il) for any j, the iterations Q apphed with initial value 
/■ = PjQ^-^Xjf will converge anyway to a function /^ with the property 

We shall now formulate in the following a reconstruction algorithm whose convergence is in 
fact a corollary of Theorem 12.21 and it is based on sampling at an almost regular (i.e., it can be 
understood as a perturbation of a uniform/regular) set of nodes. 
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Corollary 2.3. If V = {VjjjgN is a multiresolution analysis in LP{Cl) for p £ (l,oo) with prop- 
erties MRA1)-MRA3), then for any j e N there exists a j' > j large enough such that for any 

sampling set of the type S^v = {Cj',fc}fc=o '^'^2/ function fj G Vj{Q) can be exactly recovered from 

its samples {/(0',fc)}/c=o ^V ^^^ following iterative algorithm 

(8) /i"+'^=P.^*,,.H,.,(/,-/i"') + /i"\ n>l, /f =P,5*^,^H,,/,. 

In fact one has 



An) 



(9) f, :- lim f">, 

n — ^oo ■' 



where the convergence is uniform on f2, and then it is valid in L'^{Q) for all q G [1, cxd]. Moreover, 
for all f e LP{n), let fj = Pjf, one has 

(10) /= lim lim /|"\ 



J — ^oo n — ^oo 



— ^oo -' 



I JV^' 



being f!j"' defined by the iterative scheme (jSj) by means of suitable sampling values {fj{£,j\k)}k£o- 
Before working out the proof of Thereon 12. 21 we need to show the following technical lemma. 



Lemma 2.4. Under the notations and the assumptions of Theorem, \2.'A for any f G Vj{^), j £ N, 
and any d > we define the oscillation function by 

(11) oscif){x)= sup \f{x)-fiy)\. 

* yeBs{x)nn 

Then, one has 

(12) II osc(/)|U<Q<52^W2+i)|f]| 11/11^. 



Proof. Let us assume / e Vj{^). Therefore, f{x) = ^k=o^k{S)4'3,k{x) and, if y G Bs{x), then 



\f{^) - f{y)\ < J2 Mfm^^x) ~ 0,, fe(2/)| 

fc=0 

< \\Vct>j,k\\o.\x^y\J^\ak{f)\<C62^^''^'+'^Y.\ak{f)\ 

fc=0 fc=0 

< (1 + 7V,)i/2^2^W2+i) J2 |afe(/)P < CBil + iV,)i/252J('^/2+i) 

< C,52^'W2+i)|f] 



/ N- \^/2 

The estimation [J2k=o l'^fc(/)l^) — ^ll/lb for some _B > is valid because $j is a basis for 

Vj{n). Here we have denoted Cj = CB{1 + Njf/'^. Therefore, (O is valid. D 

Now we have all the necessary ingredients to prove Theorem 12. 21 

Proof. We claim that for Aj > large enough and for any A j -dense sampling set Xj := X{Aj) 
one has 

(13) II (/ - PjQ^„x,)fAoo < ^||/,||oo, for all f, G V.iil), 

for some < ?7 < 1. This would imply that PjQ^ .x is an invertible operator on Vj{fl) with 
inverse 

oo 

(14) {PjQ^„x,)-' = J2(^-p,Q^,.x,r, 

SO that 

oo 

(15) /, = {P,Q^,^x,)-'P,Q^„xJ, = E(^ - P,Q^,.x,rPiQ^.,xJj, 

n=0 
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with uniform convergence. Of course, it is immediate to see that |(SJ| and 10 are just a reformulation 
of 113). So, let us show (O and ifTH) : 

\\{I - PjQ^„X,)fj\\oo = ||P,(/-Q*,,X,)/j||oo 

< q\\il~ Q^„X,)f Aoo 

\l=0 J 1=0 

Mi 



< c;'iiEl/.(-)-/.(^.,^)lV'^ 



e=o 

Mi 



< C'^\\J2\o^scf,{x,j)\iM 



e=o ' 

< Cj||0SC/,||o,. 

By applying Lemma 12.41 one finally has 

\\{i-p,Q^,.x,)f, Woo < c,s,2^<^'/^+'^M\fAoo. 

Therefore, if Aj = SJ > is large enough, then one immediately has (|14|l . The rest of the 
Theorem is shown as a consequence of MRA3). Moreover, if / e C'(il), but not necessarily 
/ 6 Vj{fl) for any j, then P-jQq,j,Xjf £ Vji^) and 

p,Q^„xj - (p,g*,,x,(p,g*„x,)-^)p,g*,,x,/ 

Let us denote /(°°) := (Er=o(^ " ■P7g*,,x,)"/'jQ*,,x,/) £ Vj{n). Then one has 

(16) i^,(Q*„x,/(°°^-Q*„x,/) = 0. 

This proves the previous Remark. D 

Similarly Corollarv l2.3l is shown: 

Proof. Since we have assumed 2~^~ 'Ylik=o't'o',^ = ^ '-''^^ *^an write V'fe := 2~^~(j>j'^k and Xj',fc := 
^j',A; G supp(0j',fc), and one immediately see that Q^'.^x' = S<s>.,,s.,- Therefore an application of 
Theorem 12.21 would conclude the proof. Let us anyway observe more explicitly that for j' > j 
large enough, under the hypothesis MRA2) it is 

only if X G supp((j)ji ^k)- That means that 

Ifji^) - fji(j',k)\(t)fA^) < I osc fj{Cj',k)\(t)fAx), 

Cm2~i 

pointwise. Therefore, following the proof of Theorem 12 . 21 one obtains 

\\{I~P,Q^„x,)fA\oo < C,m2-^'2^W2+i)|fi|||/,||^, 
so that Cjm2~-'"2^('^/2+i) \n\ < 1 for j' > j large enough. D 



MASSIMO FORNASIER AND LAURA GORI 



2.1. On the multivariate interpolation problem on domains. An interesting interpretation 
of Theorem 12. 21 is the foUowing: 

If, for a sequence of data {y^j^lg, there exists / G Vj{n) such that f{xi) — yi and {xe}f^Q 
are dense enough, then / is unique. In fact if there were two functions f,g^Vj (fl) such that 
f{xi) = ye = g{xe) for aU £, but f ^ g, one would have Q*j,Xj/ = Q^j.x^g and the foUowing 
absurd consequence 

/ = {PiQ^,,x,)-'PjQ^,,xJ = {PjQ^,^x,)-'PjQ^„x,g = g. 

In other words, if the interpolation problem of the data {{xe,ye)}fLQ is solvable in Vj{il) then 
it is uniquely solvable. Unfortunately there are not simple conditions on {{xe,ye)}fLQ so that 
the interpolation problem is solvable in Vj{fl). However, as a consequence of (IIGII . one has the 
following guasi- interpolation argument: 

For any A-dense set X of sampling nodes in Q it is always possible to construct a piecewise 
constant interpolation operator 

M 

(17) ^x/:-5]/Mxo,, for all / e C(f7), 

where 

ne := {x e n : \xe - x\ < \x, - x\, Vi 7^ £}, £ = 0, ...,Af 
defines a Voronoi decomposition of il. In particular, one has Vx-f{xe) = f{xe), for all i = 0, ..., M. 
It is not difficult to see that substituting Q* ,x with Vx- in Theorem 12. 21 the result will be valid 
again. Moreover, if / = J^e^oTJ^Xni then, formally, it is Vxf{x) = f{x) pointwise. Therefore one 
has 

= Pj(V3c,/(°^) - VxJ) = P,{VxJ^°"^ - /), 
uniformly. 

3. Numerical implementation, examples, and results 
In this section we want to illustrate a rather efficient numerical implementation of the scheme 



©• 




Figure 1. Cubic B-spline basis on the interval [0,4] 



Without loss of generality we assume to work in the space Vb(^)- First of all it is important to 
discuss how to construct a possible projector Pq. A natural choice can be in fact given by 



(18) 



^0/ 



Nq 

E 

/c=0 



(/, 0O.fc)0O,; 



where $0 '■— {^o.fcjfei^o "^ L'^i^) is a biorthogonal dual basis for $0. In particular Pq would 
result as a projection from L'^{^ onto Vo(^)i of course bounded from C(SX) into Vb(r2) both being 
endowed with the sup-norm. 

Therefore the problem of constructing a good projector is shifted now to the construction of 
suitable dual bases. For example the canonical dual biorthogonal basis in Vq{^ is computed as 
follows: 
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Figure 2. Canonical dual basis on the interval fl = [0,4]. The support of the 
functions coincides with the domain ft. 




Figure 3. Sampling of a function in the cubic B-spline space Vo[0,4] = 
span{i?o,3(a; + 3 — i) : i = 0, .., 6}. Observe that the number of sampling nodes 
coincides with the dimension of the space Vq. However the nodes are strongly 
non-uniformly distributed with an almost coincidence of two of them. 



Consider the Gramian matrix G(<I>o) = {{(po.h, 4>a.k))h,k=a No- Then the rows (or the columns) 

of G(<i>o)~^ are the coordinates of the elements of the dual basis $o with respect to ^q, i.e., 



(19) 



$o = G($o)"'$o 



where here, with an abuse of notation, we have considered both $0 and $0 ^s column vectors. 
Unfortunately the canonical dual basis can be formed by functions with support on the whole 
domain Q. In order to obtain banded matrices and a more efficient reconstruction scheme, it will 
be more pleasant to deal with biorthogonal dual functions with support strictly contained in il. 
The construction of such duals will be discussed in Section 6. 




10 20 30 40 



Figure 4. Successive iterations of the algorithm for the reconstruction of the 
function from its samples as illustrated in Figure 3. One can see the slow conver- 
gence of the approximant to the function to be restored. 
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Now that we have a recipe how to construct Pq as in (|18|l . we should discuss how to implement 
numerically the general scheme ^. Recalling from (|17|l that 

Mo 

Q*o,Xo/ = X! /(^o,£)V'o,£, for all / e C{n), 

£=0 

one immediately has 

No Mo 

No /Mo \ 



X! ^fi^oi){^o,e,(i>o.k) 4>o;i 



fe=o \e=o J 

This suggests to define the notations P := [/(a;o,o), ••■, /(a:o,Mo)]^, M^0 := (("00/, ^o,fc))fc=o,...,Aro;£=o,...,Mo> 
^o" := ((/'o,fc(a;o/))fco,...,j\/o;fe=o,...,Afo: and *o'' := (0o,fe(T«))i=ieZ'i,rien;fe=o,...,JVo for t > small, 
and one can write 

(20) PQ,(P):=*o^M,^^P, 

(21) PQe(f^) := *o^M^^P. 

If /o G Vo(r2) is the function that we want to reconstruct from its vector of samples fg — 
[/o(a^o,Oi •■•J /o(2;o,A/o)]^i one can use the following algorithm to compute/to approximate the vector 
fo = [/o(''"*) '■ i £ ^'^j Ti G ri]"^ of its sampling on a r~^-dense regular grid. The convergence of the 
following discrete scheme is of course ensured by the uniform and pointwise convergence of the 
original iterative algorithm (0). 

Algorithm 1. RESTORE[n,f^] ^ f^: 

f- := PQ,(fo^); 
f-=PQ,(f^); 
i = 0: 
While i < n do 

i -.^i + l 

f=:=r + PQ,(f^-f^); 

f^:=P + PQ,(f^-P); 
od 

re cc 

Iq .— r . 

The discrete PQ-procedure is implemented as in H20|l and H21|l by suitable matrix-matrix and 
matrix-vector multiplications. If the functions defining the matrices are assumed compactly sup- 
ported (with support smaller than ft), then such matrices are banded. Therefore, in such a case, 
Algorithm^can be really implemented as a very fast reconstruction procedure. Unfortunately the 
canonical dual ^o,fc G ^(^) is not in general locally supported and different biorthogonal duals to 
compute a projector Pq should be considered. Moreover, since the convergence is monotone with 
a rate of decay proportional to A~", one has to expect of course faster results whenever A-dense 
sets X of sampling nodes are considered with A > larger and larger. Figures 3-11 show uni- 
and bi-variate examples of applications of Algorithm ^ in cases where the sampling sets are highly 
nonuniform and not much dense, and, as a consequence, with a relatively slow convergence. 

3.1. Computational cost. Each iteration of the algorithm depends on certain matrix-vector 
multiplications where the involved matrices 



PQ,(P) := ^o'^M^^P, 
PQ,(P) := *o^M^^f^ 
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Figure 5. Maximal reconstruction error in logarithmic scale for the first 50000 iterations. 




Figure 6. Pointwise reconstruction error in logarithmic scale at the 50000 iteration. 

can be both assumed to be banded (typically with bandwidth independent of the level j!) as soon 
as one uses suitable biorthogonal dual bases to define the projectors Pj. Moreover, the rate of 
convergence depends on the constant rj = r?(A) < 1 as in formula ifT^ so that 

||/-/(")||oo<rr||/||oo. 

Therefore, for a fixed density A, for achieving a prescribed fixed tolerance e > the computational 
cost can be assumed 0(#supp(f'^)), where P is the output. This means that the algorithm has 
optimal complexity. 





Figure 7. The original surface belonging to Vb([0,3]^) = span{_Bo.2(a; + 2 — 
i)BQ^2{y + 2 — j) : i, j = 0, .., 4} to be reconstructed from its samples is shown on 
the left. On the right the result of the reconstruction is shown after 6000 iterations. 



4. Multivariate iterative numerical integration 

Under the assumptions and notations considered so far, let fj G Vj{fl). The iterative formula 
(|S|) in Theorem 12 . 21 converges uniformly to fj and this implies also that 



(22) 



fj{x)dx = lim 



4"+'^(x)dx. 
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Figure 8. 100=4 dim(Vo) sparse random sampling points are chosen. 









Figure 9. Some iterations of the reconstruction scheme. 




D :oi>De laaa 3«»a taaoa waaji »d« mm 



Figure 10. The || • ||oo-error is shown in log- scale for successive iterations. 



On the other hand 



f]''^'\x)dx 



f(")/ 



PjQ^,.x,{fj-fr'){x)dx+ / f/"{x)dx 
SI Jn 



k^O £^0 



-jMyf|-(frn+//r(-)^-' 
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Figure 11. The final local pointwise error in log-scale is shown. 



where M„ 



i{'ipj,e,(f'j,k))k=0,...,Ni-i=0,. 



,Mi 



[fjiXj,o),--J]{xj.M, 



r, (fj"V 



and ojj = [/q (j)jfi{x)dx, ..., J^ 4'j,Nj {x)dx]'^ . Therefore, similarly to Algorithm^ one can formulate 
the following iterative integration procedure. 

Algorithm 2. INTEGRATE[n, f?] -^ ly. 

int := ujjM^ -ff ; 

i = 0; 

While i < n do 

I := i + 1 



[/f(x,,o),...,/i"^(x,,M,)]^, 



int := int 



- OJ, 



^M;^(fi^-f^)5 



P + PQi(ff-P); 



od 
h : 



int. 



Here PQ^ is defined analogously to H2Q(I . 



the procedure INTEGRATE has the 



Proposition 4.1. Under the assumptions of Theorem\ 
following properties: 

i) It is linear, i.e., 
(23) INTEGRATE [n, Af^ + /ig^] = A • INTEGRATE [n, f*] + a* • INTEGRATE [n, g^ 



for all X, fi £ C 
ii) for all fj G Vj{Q) it is 



(24) 



/j(x)dx==INTEGRATE[oo,fJ^] := lim INTEGRATE[n,fJ*; 



Proof. To show i), it is sufficient to observe that all the operations executed in the procedure 
INTEGRATE are linear. The proof of ii) is a straightforward apphcation of if^ . D 

REMARK: The previous Proposition u) means that integration formula INTEGRATE[oo,f^] 
is exact on fj e Vj{n). Of course, one never can compute INTEGRATE[c>o, f?] but only 
INTEGRATE[n, f?] for n large enough to ensure the desired accuracy. The accuracy is achieved 
without increasing the density of the sampling points f but just going further in the iterative 
integration process, see Table 1. 

Corollary 4.2. // Vxi^) C Va{fl), where Vxi^) denotes the set of algebraic polynomials of 
degree K Cz N at most, then INTEGRATE[oo, p**] is an exact iterative integration formula for 
all polynomials p £ Vk{^)- 
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Iter. 


Approx. integral 


Error. 





1.07874 


0.131095 


50 


0.97298 


0.0253353 


100 


0.953824 


0.00617898 


300 


0.947611 


0.0000340001 


2000 


0.947643 


1.35891x10-6 



Table 1. Application of the iterative integration. Integral value « 0.9476446652462314. 



5. Error analysis 

As we have shown in Section 2, for a continuous function / ^ Vj for all j > 0, Algo- 
rithm 1 applied on A^ > dense samples of / computes a function / with the property 
PjiQ^j.Xjfj — Q^j,Xjf) = 0. This is an indirect information on the error that we have in the 

approximation / w / . In this section we address the direct estimation of the error ||/ — /, ||2, 
under the assumption that / has some regularity. In particular we will assume that / £ if (fi), 
for s > 0, where H^iVt) denotes the Sobolev space of order s. Moreover, we also estimate the error 
that one produces by approximating the integral of / by INTEGRATE[cxd,P]. The key tool is 
the characterization of H''{il) by suitable biorthogonal wavelets. In particular we assume in the 
following that 



Wl) there exists a wavelet biorthogonal basis {'^j.k}j>o,keJ associated to the MRA {Vj (il)}j>o, 
so that any function / G i^(i7) can be written 

(25) /- Y. if'hk)'l^j,k+ Y. {f,^j,k)^j,k^Pjif) + f^; 

k=0,...,Nj j>J,keJj 



W2) / eiJ*(n), s > 0, if andonly if 
(26) 



1/2 



1/2 



ci\\f\\HHn)<\\Pjfh+\ J2 22^-'ia*.,fc)P <c2\\f\\HHn), 

\]>J.keJ, J 

for ci, C2 > independent of /. This characterization of H^(ri) ensures in particular that 

\ 1/2 / \ 1/2 

< 2-M E 22«^|(/,§,-fc)p 
\]>J,keJj 

W3) ||*j,fe||oo < C^2^ and #J,- < Cj2*'; 
W4) /j-^ *j,fc(a;)dx = for aU j > and fc G J,. 

Proposition 5.1. Let f G C{fl)r\H''{fl) for s > 0, but not necessarily f E Vj for any j > 0. Under 

the assumptions of Theorem \2.'A one has that / := Sjf = (X]^o(^ ^ PjQ<i j ,XjY^ PjQ^i j.Xj f) G 
VjiVl) has the property 

(27) 



„(oo) 



\\f-fr'h<C\\S,\\2-^^\\f\\u^^a). 

where \\Sj\\ is the norm of the operator Sj on i^(il). 

Proof. With similar arguments as in the proofs of Lemma 12.41 and Theorem 12.21 one shows that 
the map Sj : f —^ /, is bounded with the L^-norm. Moreover this map coincides with the 
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identity on Vj. This together with the properties Wl-2 imply 



11/-/ 



(oo). 



= \\f-S,if-P,f + P,f)h 

= \\Sjif-P,f)-{f-S,{P,f))h 

= \\sAf-Pjf)-if-Pif)h 

< (1 + 11^,11)11/ -P,/||2 

< C'(l + ||5,||)2-^^|l/||^=(n). 



n 



Lemma 5.2. Let f G C{fl). Under the assumptions of Theorem \2.<^ one has that 

(28) |INTEGRATEfcx),P]| < ■i^^^J^II/lloo, 

1 — 1] 

where < rj < 1 is as in formula 1)13(1 . 

Proof. By the remark in Section 2 and by ((22|l . one has 

|INTEGRATE[oo,P]| = 



< 11^1 



oo 

Y.{i~p,Q^,,x,rPjQ^,„xj 



n=0 



< 



< 



\m 



\PjQ^,,xJ\[ 



1-77' 

\mPMQ^„x,\ 

1 — rj 



One concludes observing that HQ* ,x || = 1- 

Then one has the following error estimation result. 



n 



Proposition 5.3. Assume that ^ < s. If the iterative integration formula is executed on a Aj- 
dense set Xj , for J > and for Aj > large enough, then there exists a constant C > such 
that, for all / £ H''{Q), it is 

2(M-s)(,7+l) 



INTEGRATE[cx),P] - / f{x)dx 



< rMMiini 



1-2^ 



Proof. First of all observe that s > | and then, by the Sobolev embedding theorem, it is / G C{fl) 
and it makes sense to consider its sampling. By Proposition 14. II il that ensures the linearity of 
INTEGRATE and by Wl) 

INTEGRATE[oo, r] - / fix)dx 
Jn 

INTEGRATE[oo,Pj(f)=] + INTEGRATE[oo, V (f, *j,k)*J,k] - / Pj{f){x)dx+ V {/,*,, / *,,fc(a;)da:| 

j>j,k y^ j>j.k ■'n / 

In the last equality we could exchange the integral with the sum because the sum converges 
absolutely due to W2) and W3): 

j>J,k j>J.k j>J.k j>J 

Moreover, since INTEGRATE is exact on Vj by Proposition 14. II ii) and J^'^jk{x)dx — one 
has 



(29) 



INTEGRATE[oo, P] - / fix)dx 

n 



INTEGRATE[oo, E (f' *j,k)*j-,k] 
j>j,k 
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By Lemma 1^1 equation (f^ . and by using again W2) and W3), it is 



INTEGRATE[oo, P] - / f{x)dx 



< 



M\Pj\\ 



< c 



c 



1"^ 

'mpj 



l-r, 

\mpj\ 



j>j,k 

2(f-s)(J+i) 



H^ 



1-2^ 



n 



REMARK: In other words, the previous theorem tells that if / £ H^(i}) with s > ^ and one 
has just a Aj-dense set of sample points {fixj,i)}^Q, for Aj large enough, then its (multivari- 
ate) integral can be computed by using the procedure INTEGRATE with an error of order 
0(||P/||2'^~^'*^^-^+^^), once a sufhcient number of iterations has been done. 



6. Construction of MRA on domains and manifolds 

In this section we want to recall the construction of a rather large class of multiresolution analy- 
ses with properties MRAl-3 for which the axiomatic sampling theorems illustrated in the previous 
sections can be applied. It is also ensured that corresponding wavelet bases with properties Wl-3 
are also available. 



6.1. Composite multiresoultion analyses. The main ingredient is the construction of biorthog- 
onal refinable bases $ ' :— {0j,Oj ■ ■ ■ , 4'j.Nj } ^-nd $ ■ ' := {0j,O) • • • » 4'j,Nj } to define the multiresolu- 
tion spaces with complementary boundary conditions on D = (0, 1)"^'. In 12 such bases have been 
derived from integer shifts of B-splines On on the real line. First one considers a biorthogonal dual 
bases 9^^ ff as described in 0, and then both Qn and 6^^ jv ^^^ restricted to define bases on (0, 1). 
Those elements of Qn and 6^ fj that do intersect the boundary of (0, 1) are modified in order to 
ensure i) the polynomial reproduction with maximal degree, and ii) to match prescribed comple- 
mentary boundary conditions. Such modifications of course imply that the resulting bases are no 
more biorthogonal. However such new bases can be in turn biorthogonalized. Finally, the bases 
<& ■ ' and $ • on D are defined as tensor products of the resulting biorthogonal dual bases on (0, 1). 
In |Hl|71^1Ej the construction of spline- type multiresolution analyses on rather general manifolds 
has been proposed. Such construction is based on the decomposition of the manifold Q into sub- 
manifolds Jli smoothly parametrized by D := (0, 1)'', i.e., fl = uflifli, fli — Ki{0), i — I, ...,M, 
where fti Dftj — ^ for i ^ j, and k^ : R'' — ^ R'* , d < d' are smooth regular functions. The idea is 

to define biorthogonal multiresolution analyses V, ' (D) and V"- ' (D) on D adjusted with suitable 
complementary boundary conditions Zi and Zi as described above ^21) for each i = 1, ..., M . Each 
MRA is lifted on fJ,; by using the parametrization k^ to define Vj{p.i) := {f o k^^ : ip E V, '(□)} 
and similarly Vj(Qi). The boundary conditions are chosen to fit with a suitable bounded resolution 
of the identity 7^ = {Ajfii on n so that Vj{n) := Y,fi^ i?jV,(f7,) and Vj{n) := Y,^^ R-VjiQ^) 
will define global biorthogonal multiresolution analysis on H. with properties MRAl-3. Moreover, 
associated to such global multiresolution analysis one can also define suitable biorthogonal wavelet 
bases ^3] with properties Wl-3. 

Unfortunately the definition of the projectors Ri makes use of certain extension operators that 
show increasing norms as soon as higher smoothness is required. This makes the resulting bases on 
n rather ill conditioned. One way to compensate this drawback is to start at the very beginning 
with refinable functions that are better conditioned than B-splines. 
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Figure 12. GP functions are illustrated for n = 2 and n = 3 and for different 
choices oih> n — \. 

6.2. GP bases. A more general class of positive compactly supported refinable functions on the 
real line has been characterized in |52] as the solutions of the refinement equations 

n+l 



(3U) 








'^(•^) = Z^«fc' ^(2a 










fe=0 


where 










(31) 


4-'^) 


^2-^ 


A fc ) 


+4(2-"-i)(^:;)' 



?i>2,/i>n-l,A; = 0,...,n+l. 

Let us call GP-functions the elements of this class. In particular B-splines are GP-functions for 
h = n. 




Figure 13. Normalized GP functions adapted to the interval (0, 1) [^ Section 
4] for n = 5 and h — 10. 

This generalization is worth for the following reason: It has been pointed out in j^Jj (see also 
|2()ll23j) that they are better conditioned as bases O^"'''-' := {ip{x — k) : fc G Z} for increasing values 
of h. This is due to the remarkable properties to have a smaller and smaller essential support, see 
Figure 12, while preserving the degree of smoothness, as soon as h increases, see Table 2. 



h 


3 


4 


5 


6 


7 


8 


20 


K2 


46.40 


26.06 


20.53 


18.68 


18.00 


17.73 


17.49 



Table 2. The spectral condition number K2 of the Gramian matrix for GP 
functions adapted to the interval (0, 1) for n — 3 and increasing values of h. 
In particular, in |21| it is shown how such bases can be adapted to the interval together with 
certain biorthogonal duals as derived in [^ although no complementary boundary conditions have 
been considered yet. We postpone this second adaptation to a successive contribution. 
Let us shortly show instead a different method to produce biorthogonal duals for GP-bases as a 
generalization of the results in [S|. As we have already stated in Section 3, the construction of a 
compactly supported biorthogonal dual is fundamental for the definition of suitable projectors Pj 
raising banded matrices in Algorithms 1 and 2. 
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Associated to a mask a = {a^ : k = ni,...,n2} for the refinement equation, one defines the 
symbol 

(32) p{z) := Y. «''^''- 

k—Tii 

Equivalently one defines the trigonometric polynomial 

(33) m(0 ■■= \p{e~'^)- 

A refinable function ip with mask a'-"'''' — {flj"' : k = ni,...,n2} defines a biorthogonal dual 
basis G)("''') := {(^(x — k) : k £ Z}, i.e., J^if{x)ip{x — k)dx — Sk,o, only if its associated symbol 
ffii-n-M) Y^Qjg f^ijQ property 

(34) TO("''')(Om("^''HO + "^'"'''^(C + 7r)m(«^'')(e + 7r) = 1. 

Therefore, one can try to find a suitable trigonometric polynomial m!-'"''^^ satisfying (|34|l . One has 
first the following result. 

Lemma 6.1. For n > 2 and for h > n — 1 it is 

Jn,h)(C\ - ^-i^i o^^ f i 1 1 /,'o'i-n+lV9'i-«+l 



(35) m^"'''J(0 = e-'T-'^cosl - 1 l/(2'*-"+^)(2'^-"+^ - 1 + cos(C)). 

Proof. In 22, Formula (3.9)] it is shown that 

/"^'')(z) = 2-''(l + z)"-i(z2 + (2'^-"+2 _ 2)z + 1). 
A direct computation shows that ni^"^'''' (^) :— ip*^"'''^(e^*^) is as in (|35(l . D 



Therefore, to'"''''(^) = e ' 2 ? cos ( | J r(cos(a;)), where A'^ G N and r is a certain polynomial. 
One can look for solutions of (|34|) of the same form, i.e., to'^"''')(^) = e^*T~^ cos ( | ) f(cos(a;)), 



2^ 

where 2^ := A^ + A^ is even and f is a suitable polynomial. Substituting these expressions in H34|l 
one obtains 

(36) *^°^ ( 2 ) ''('=os(^))^(cos(0) + sm f - j r(- cos(C))r(- cos(O) = 1- 

This equation is exactly the same as |S1 Formula (6.9)], and it is known to have solutions f 
constrained by 

r(cos(0)r(cos(C)) - E (^ " • ^ ') ^^^ (|) ^ + «"^ (f ) ^(cos(C)), 

where R is an odd polynomial. In other words one looks for a polynomial f{x) — a^ + aix + ... + 
amZ™ such that 

-1 + A/i-a; 



s[ao,...,am\{x) ■.= r{x){ao + aix + ... + arax"^) -'^i . jf 

is a polynomial divisible by (— ^)^ with an odd polynomial as the quotient. Since r{x) — 
l/(2''~"+^)(2''~"+^ — 1 + x) is a polynomial of degree 1 and the quotient R has at least de- 
gree 1, it is clear that f must be a polynomial of degree m at least i. Let us assume then m > £. 
Therefore one reduces the problem to the solution of the following system of equations 

mod (s[ao,...,a„](x), (i^) j = 0, 

s[ao,...,am]ix)/ {^Y + s[ao,...,a„](-a;)/(i;^^) =0, 
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where p{x)/q{x) here indicates the quotient of the division of the polynomial p by the polynomial q. 
Since the remainder mod (s[ao, ■■■,am]{x), (^3^)^) is a polynomial of degree £—1, the first equa- 

tion will impose £ linear conditions on the aq, ..., a^n variables. The quotient s[ao, ..., a,n]{x)/ (— ^) 
is instead a polynomial of degree m + I ^ £. Therefore the second equation imposes further 
L "'"^2^~^ J + 1 equations. Altogether one generates [ '"^^^~^ J + i + 1 conditions. In particular for 
m = £, m = £ + 1, and m = £ + 2 one has that the number of conditions is exactly m + 1, as 
many as the number of unknowns uq, ..., am- For m > £ -\- .i the number of conditions is strictly 
less than 771 + 1. This allows to impose further conditions to determine uniquely the solution. 
Once determined a trigonometric polynomial ifv-'^'^ by solving the linear equations, one should 
check that the equivalent conditions 8, Theorem 4.3 Cl-3] are satisfied in order to conclude the 
successful construction of an admissible biorthogonal dual. 

Example 1. Let us consider n — i and h = n + I = i. By Lemma [6.11 it is 



Since in this case iV = n — 1 = 2, we choose N = 2 and 2£ = N + N = 4. Therefore one has to 
look for a polynomial f{x) = ao + aix + a2X^ of degree at least 2, such that 

r ■] / \ ^ ' "^ / 2\ \ ^ I ^ ' ^\ /t X 

s[aQ, ai,a2\(x) = — - — (ao + aix + a2X ) ~ 2_^ 
is divisible by ("3^) , i.e., 

mod s[aQ,ai,a2]{x), ( ^—^ ] I = -2+3/4ao-l/4ai-5/4a2+(l+l/4ao+5/4ai+9/4a2)a; = 0, 





and its remainder is an odd polynomial, i.e., 

fl~x\'^ /l-(-a;)V 
s[aQ,ai,a2]{x)l [——\ + s[ao,ai,a2\{-x)/ [ j = 2a2 + lOas = 

These conditions are equivalent to the solution of the following linear problem: 

3/4 -1/4 -5/4 \ / ao \ / 2 
1/4 5/4 
2 

The solution is (00,01,02)'^ = (8/3,-25/12,5/12)'^ 

to(3'4)(^) ^ g-2^« cos(C/2)2(8/3 - 25/12 cos(C) + 5/12cos(e)'). 

The corresponding mask is therefore given by 

x(3,4) _ r„(3,4) _ 5 (3^4) _ 5 (3,4) _ 43 (3,4) _ H (3,4) _ 43 (3,4) _ 5 (3,4) _ 5 

that defines again a symmetric filter, see also [HI. A direct numerical computation shows that 
j^(3,4) gp defined verify together with 7tt,(3>4) ^j^g conditions 8, Theorem 4.3 Cl-3] so that 

5 

^{x) = ^ a^''*<^(2a; - fc), 

is an L^(M) function and ((p, (f{- — k)) = Sk,o for all fc G Z. 

As already mentioned above, the composite biorthogonal bases as derived on manifolds in 
|13j might exhibit limited smoothness or high condition numbers. As proposed here, the use of 
better local bases certainly improves this drawback. However an alternative way to overcome this 
difficulty is to use frames instead of biorthogonal bases, see for example |11II25| . Such frames are 
constructed by Overlapping Domain Decompositions where the patches are again smooth images 
of D. The fact that one does not need to implement interfaces through patches to preserve global 
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smoothness reduces the ih conditioning of the global system. The use of frames do not affect, e.g., 
the core of the proofs of Theorem 12 . 21 and Proposition 15. II 

6.3. Implementing GP functions for the sampling problem. As stated in the previous 
subsection GP-bases 6("'''' adapted to the interval have better condition numbers for increasing 
values oi h > n — 1. Thus one is tempted to affirm that for any / e V^("''*) := span 0("''^) for 
all h > n ~ 1, f will be reconstructed with increasing rate of convergence from a fixed A-dense 
sampling set, for increasing values of h. In particular, it has been shown |55| that Vn-s C l/("'^) 
for all h > n — 1 and, therefore, one can test this claim on polynomials. Surprisingly the claim is 
false. Numerical experiments for the reconstruction of the polynomial f{x) — {x — i)^ — 0.1 on 



show 



(0, 1) from a nonuniform sampling based on Algorithm 1 implemented by using V ' , j — 5, 
that the asymptotic rate of convergence decreases for increasing values of h. However, this also 




Figure 14. Nonuniform sampling of f{x) = {x - i)^ - 0.1. 

means that by using V ' , j = 5, and n — 1 < ft, < n one can anyway obtain better performances 

in the reconstruction than implementing Algorithm 1 by using just the B-spline space V ' , see 
Figure 15. 
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Figure 15. The solid line represents the uniform logarithmic error in the re- 



construction of f{x) 



l\2 
2' 



0.1 by using the B-spline space Y- 



(5,5) 



The 



dashed lines represent the uniform logarithmic error in the reconstruction of 
/(x) = {x- \f - 0.1 by means of vf'^'^ for h G {4.1,4.3,4.5}. One can see 
that for h = 4.1, one achieves the machine precision error with a 20% less of the 
number of iterations than for the choice h — n — h^ corresponding to B-splines. 



6.4. Changing parameters in quasi-interpolation. For a given sampling set {(x^, yi)\\LQ and 
a fixed space V,, it is not ensured that there exists /j G Vj such that fj{xe) = yg for i? = 0, ...,M. 
Nevertheless, as soon as the set is dense enough Algorithm 1 will converge anyway and the resulting 
function / will depend on several different parameters: The choice of the quasi-interpolation 
operator Q<s/^x and the choice of the space V}. Therefore one can tune the choice of Qip^x and Vj 
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in order to optimize certain geometrical properties of the resulting / with respect to the given 
data set {{xe,ye)}e'Lo: ^^^ ^^^ example Figure 16. 






Figure 16. Segmentation of the brain phantom. The curves are computed with 
Algorithm 1 from the nonuniform sampling indicated by the white pixels. Dif- 
ferent curves are generated by changing the GP space V and the quasi- 
interpolation operator Q^Sf.x- 
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